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Abstract. Time-dependent nonequilibrium Green’s functions (TDNEGF) are shown 
to provide a flexible, effective tool for the description of quantum mechanical single 
particle scattering on a spatially localized, time-dependent potential. Focusing on 
numerical methods, arbitrary space and time dependence of the potential can be 
treated, provided it is zero before an initial time instant. In this case, appropriate 
version of the Dyson and Keldysh equations lead to a transparent description with 
clear physical interpretation. The interaction of a short laser pulse and an electron 
propagating initially in free space is discussed as an example. 
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1. Introduction 


Time-dependent scattering problems appear e.g. in atomic, molecular and solid state 
physics, but they are of relevance also in nanoscale electronic devices driven by an 
alternating or pulse-like bias. In spite of the difference between these physical systems, 
the mathematical models that can describe these phenomena are similar. In the 
current paper we focus on methods based on nonequlibrium Green’s functions, which 
- besides providing an adequate description of many particle quantum systems - 
have important applications also in transport-related problems in solids |T|. Ballistic 
(coherent) transport in the single electron approximation (which is often used for 
nanometer-sized samples) is equivalent to single particle quantum mechanical scattering 
problems. That is, methods developed for the description of transport processes in solid 
state systems have a more general scope in scattering theory. 

The formalism based on time-dependent, nonequilibrium Green’s functions 
(TDNEGFs) has been initiated in the early 1960s 12,3 but has many fundamental 
aspects that are discussed also in recent textbooks on field theory j4j[5]. This approach 
has been successfully applied to transport phenomena mainly in nanoscale solid state 
devices and static potentials (for a recent summary, see (6|). TDNGF-based description 
of time-dependent transport in mesoscopic systems was considered already in 1993 |7|, 
but numerically effective methods have been developed more recently. The details of 
these methods can be found in reference [§], where the connection between TDNEGF 
and other approaches [9,101 is also discussed. 

In the current paper we investigate how the formalism based on TDNEGFs can 
be used to describe scattering phenomena. For the sake of simplicity, we consider one 
spatial dimension, which is useful for the clear interpretation of the final formulas, but 
does not mean a necessary restriction. (Generalization to two and three dimensions 
can be done without essential difficulties.) We consider a monoenergetic plane wave 
that propagates initially in free space, i.e., the potential is zero. Then, at t = 0, a 
localized potential emerges, causing e.g. time-dependent reflection phenomena. In order 
to keep the generality of the treatment, a spatial grid is applied, i.e., the position 
variable is discretized. Note that this is the only numerical approximation in the model 
(allowing calculations for potentials with arbitrary space and time dependence), all other 
calculations leading to the dynamical equations are analytic. 

Time-dependent theory of nonequilibrium Green’s functions rely on the Dyson 
and Keldysh equations as central concepts, and usually the main goal is to determine 
the so-called lesser Green’s function from which direct physical consequences can be 
drawn. The solution of the equations resulting directly from the TDNEGF theory mean, 
however, generally an extremely complex task, the numerical costs of such calculations 
are tremendous. In our case, as it is going to be shown, the combination of the Dyson and 
Keldysh equations leads to a time evolution that can be calculated efficiently, and have a 
clear physical interpretation as well. Namely, the dynamics has to be solved only in the 
spatial region where the potential is nonzero, in this domain the usual time-dependent 




Quantum mechanical scattering on time-dependent potentials: TDNEGF method 3 


Schrodinger equation is valid, and TDNEGF theory takes the boundary conditions 
into account by allowing the wave function to leave the ’’region of interest” without 
disturbance. These ’’perfectly transparent” boundary conditions turn the Schrodinger 
equation into an integro-differential equation, but the relevant integrals involve only the 
boundary points of the finite region in which the dynamics is to be calculated. Note 
that these results are analogous to those derived in reference |8| in a different way for 
two dimensional ballistic transport (and applied for the description of ’’flying qubits” in 
coupled quantum wires), but here the focus is on general scattering phenomena and the 
interpretation of the dynamical equations. Additionally, as an application, in the last 
section we calculate how a short laser pulse interacting with a monoenergetic electron 
beam can create electron density fluctuations. 

2. Statement of the problem 

Let us consider the problem shown schematically in figure [TJ For t < 0, we have a 
monoenergetic quantum mechanical plane wave propagating in free space 

$(x,t <0) =e i ^ E)x - Et \ (1) 

Here, we set ft = 1. (This convention will be used throughout the paper.) According to 
the Schrodinger equation, 

k(E) = V2mE, (2) 

where m denotes the mass of the particle, and the direction of propagation has been 
chosen to correspond to figure [lj Then, we assume that a localized potential, the support 
of which is the interval [0, L\, is turned on at t — 0. That is, 

U (. x , t) — 0, if t < 0, or x ^ [0, L\, (3) 

otherwise the dependence of U on space and time can be practically arbitrary. 

In the following we present a solution to the problem described above. Clearly, it is 
not the only possible approach, but the TDNEGF-based method treats the scattering¬ 
like boundary conditions in a specific, exact way. In ID, less sophisticated brute force 
numerical methods can solve the problem, but they are difficult to transfer to two or 
more dimensions (which is quite straightforward using TDNEGF [8]). 

3. Green’s functions and discretization 

In single particle quantum mechanics, we can call a function with two spacetime 
arguments Green’s function, if it is the ’’inverse of the Schrodinger operator,” or, in 
other words, it is the fundamental solution of the time-dependent Schrodinger equation. 
That is, 

- H(x u h)] G(x i, £i, x 2 , f 2 ) = S(ti- t 2 )S(x 1 - x 2 ), (4) 
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a): t<0 



Region III 





Figure 1 . (Color online) Schematic representation of the problem we consider. 


which, however, does not define G(xi, t\, x 2 , t 2 ) uniquely. Probably the most important 
member of this class of Green’s function is the retarded propagator, which satisfies the 
additional condition 

G R (x 1 ,ti,x 2 ,t 2 ) — 0 for t 1 <t 2 . (5) 

As it can be shown [4], G R indeed propagates the wave function of the system forward in 
time, which follows from the fact that it essentially equals to the position basis matrix 
elements of the time evolution operator: 

G R (x 1 ,ti,x 2 ,t 2 ) = : ~-i(x 1 \U(ti,t 2 )\x 2 ) for t l <t 2 . (6) 

As we shall see later, it is useful to embed these functions in a wider context, and 
consider the one particle subspace of a many particle Fock-space and express G R in 
terms of the expectation value of field operator products: 

G fl (l, 2) = G r Oi ,ti, x 2 , t 2 ) = -iG(t 1 - t 2 ) ( [^(xi, ti), t 2 )] ) • (7) 

The commutator above generally corresponds to the case of fermions, but we do not need 
to distinguish bosons and fermions when considering a single particle. The expectation 
value has to be calculated in the Heisenberg picture (’’static”) state of the system, which 
is generally represented by a statistical operator, but in our case it simply corresponds 
to the single particle state ([!]). The Heaviside function 0 ensures that the retarded wave 
function disappears for t\ <t 2 . 

The theory of (many particle) nonequilibrium Green’s function also uses the lesser 
Green’s function (two point correlation function), which can be written as 

G K { 1 , 2 ) = G < (xi,ti,x 2 ,t 2 ) = i <'h t (x 2 ,t 2 )'b(^i,n)) • 


( 8 ) 
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These functions can be used to calculate physical quantities, e.g., the local particle 
density reads 

n(x, t) = — iG < (x, t, x, t). (9) 

[Note that the density above is not normalized in the sense that its integral over the 
whole infinite spatial domain is not unity - for the boundary conditions that we consider 
this kind of normalization is not even possible. Instead, we chose the convention 
n(x,t) = 1 for the monoenergetic plane wave given by equation ([TJ).] 

In the following we are going to show how many particle TDNEGF dynamical 
equations (that describe the time evolution of G R and G K ) give rise to physically 
transparent equations in the one-particle subspace leading to an appropriate method 
to treat time-dependent scattering problems. 


The equations of motion for G R and G < can be obtained using a method developed 
by Keldysh and others |2,3|. Since our aim is to focus on the general numerical solution 
of these equations, we recall their discretized version here, i.e., we consider a spatial grid. 
Let the stepsize be denoted by a, i.e., we have Xj = ja as grid points (j — 0,..., M for 
region II in figure [I]). Then the Green’s functions turn into matrices 

G R j(t 1 ,t 2 ) = G R (x i ,t 1 ,x j ,t 2 ), ( 10 ) 

Gf j {t l ,t 2 ) =Gf j {x i ,t 1 ,x j ,t 2 ), ( 11 ) 

and the input plane wave ([TJ also has to be replaced by its discretized version 

$j(t < 0) = $( Xj , t < 0) = e i[k{E)ja ~ Et] . (12) 

For the sake of simplicity, the Laplacian appearing in the one-particle Hamiltonian 

H = H° + U(x, t) = —— A + U(x, t), (13) 

2m 

is usually represented by a simple three-point finite difference operator, leading to 

U(xi,t) + -K, if i — j = 0, 


= Hfj + Uij(t) = 


1 


2 ma 2 ’ 


0 


if \i — j\ = 1, 
otherwise. 


(14) 


The discrete version of the ’’dispersion relation” corresponding to this Hamiltonian reads 
1 


E = 


mar 


-(1 — cos ka ), 


(15) 

thus k(E) in equation ( Jl2j is given by 

k(E) = arccos (l — Ema 2 ) . (16) 

[Note that the leading term in the series expansion of equation (15) is k 2 /2 r m, in 
accordance with ([2j.] 


Despite discretization, the problem is still infinite dimensional: the grid points ar¬ 
range from — oo to oo. More specifically, the integer index j has values in the interval 
(—oo, —1] for region I, [0, M\ for region II, and [M + 1, oo) for region III. However, it 
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can be shown [i ll that it is sufficient to concentrate on the spatial domain where 
the potential is nonzero, i.e., on region II. (For constant potentials, simple matrix 
multiplication leads to this result, but it holds also in our case.) The effect of the 
semi-infinite domains I and III are taken into account by so-called self energy terms 
(memory kernels) in the Dyson equation: 

ti 

z^G*(f 1 ,f 2 ) = H(f 1 )G i? (f 1 ,f 2 ) + [ £ r (H, t)G R (t, t 2 )dt, 


dt i 


(17) 


t2 


where boldface symbols emphasize matrix properties of the objects (the products of 
which are meant to be matrix multiplication). Note that due to the retardation of the 
functions in the integrand (t\ > t, t > f 2 ), the limits of the integral can be extended to 
±oo. The nonzero elements of H R are related to the semi-infinite domains: 


£S(M') = 


d2 9io(t, t') if i = j = 1 or i = j = M, 
0 otherwise, 


(18) 


where g R (t,t') denotes the retarded Green’s matrix of a semi-infinite discretized line 
segment not being coupled to the central region (II), and d = is the off diagonal or 


’hopping’ matrix element of the Hamiltonian (14). The index 0,1 means that one needs 


the matrix element between the termination point (indexed by 0) and the only neighbor 
it has. 

In agreement with equation (|5]), the initial condition for the integro-differential 
equation ( fl7| ) is 

lim Gf k (t 1 ,t 2 ) = -i8 jk . (19) 

Cl \C2 


Since Y, R can be calculated analytically (see the next section), equation (17) together 
with (19) allows us - at least in principle - to calculate the retarded Green’s function 


(matrix) for arbitrary two spacetime points. However, according to the traditional 
treatment, one also needs to calculate the lesser Green’s functions, in order to obtain 
physically interpretable results. In other words, the Keldysh equation 

tl <2 

i-R/j. j./\ j. j./\lt 



G < ((i,i 2 )= / / G R (t u 0 [G r (« 2 , (')] **' 


( 20 ) 


—oo — oo 


has to be solved as well. The lesser self energy S < generally corresponds to scattering 
from regions I and III into the central one. However, for the boundary conditions that we 
consider (incoming monoenergetic plane wave from the left, i.e., from region I), region 
III has no influence. Similarly to the matrix S < can be expressed using the lesser 
Green’s function g < of a ’’standalone” semi-infinite line segment: 


*') = M*o d2 &t, 0- (21) 

The analytical form of the single nonzero matrix element will be given in the next 
section. 
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4 . Solvable, numerically exact dynamical equations 


In this section, first we give explicit expressions for the nonzero matrix elements 
of T, b and S < , discuss the structure of the Dyson and Keldysh equations, and 
finally introduce numerically exact dynamical equations that can be solved without 
considerable computational costs. 

The key observation for the determination of the Green’s functions of the semi¬ 
infinite line segments (region I: input, region III: output) is that since the Hamiltonian 
acting on these domains has no explicit time dependence, time translation is a symmetry 
of the problem. Practically, g R (t,t') = g R (r) and g < (t,t') = g < (-r), with r = t — t!. 
This allows us to perform a Fourier transformation with respect to the time difference 
r, the conjugate variable of which is the energy (h = 1). The retarded Green’s function 
of a semi-infinite line segment in energy domain can be calculated e.g., by eigenfunction 
expansion and contour integration jTJ|TT| . The result is 

S$,(E) = -i e“< E >“, (22) 


where a is the spatial step size and k(E) is given by equation (16). That is, 



1 

d 


1 — Ema 2 + i\J 1 — (1 — Ema 2 Y 


(23) 


The inverse Fourier transform of this expression can be calculated analytically, leading 
to: 

sSW = (24) 

where first order Bessel’s function of the first kind appear. Since in our case S < (r) = 
d 2 g< (r) describes monoenergetic scattering into region II, its time dependence is 
particularly simple, it contains a single Fourier component: 


S< (r) = id 2 N(E)e~ iET 5 kl 5 k0l (25) 

where E denotes the energy of the incoming plane wave ([Tj) . As we shall see later, the 
convention that the constant particle density that corresponds to equation ([!]) is unity, 

implies the normalization: N(E) = \J 1 — (l — ^) 2 . 


Being equipped with the explicit expressions 


and (25), we can turn to the 


analysis of the structure of equations © and ( [20| . First, let us note that the retarded 
Green’s function (in our discrete case: matrix) of region II, the dynamics of which is 
determined by the Dyson equation is general in the sense that its knowledge is 
sufficient to determine physical quantities corresponding to different possible boundary 
conditions given by S < (r). Although in the current paper we are not going to investigate 
boundary conditions other than monoenergetic input from region I [that corresponds to 
equation (25)], it is worth emphasizing that this is not the only possibility. 
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The structure of the equations to be solved in order to obtain e.g. the local particle 
density is summarized in figure [2] Using equations (25) and (21), it is easy to see that 


if we are to obtain G^(t,t), the Keldysh equation (20) factorizes: 


Gfi(t,t) = id 2 N(E) 


G R 0 (t,t')d 2 e lEt 'dt' 


(26) 


Additionally, only the 0th column of the matrix is needed, and as one can check easily 
by matrix multiplication, the Dyson equation (JTt]) is closed to these matrix elements 
(i.e., no other column than the 0th appear on the right hand side, provided this holds also 
for the left hand side). That means that the problem related to a (M + 1) x (M + 1) 
matrix reduces to the dynamics of a vector of dimension M + 1. This means that 



Figure 2. (Color online) The structure of the equations to be solved in order to obtain 
the particle density that is proportional to G5 (t, t ). 


i.e., 


for a given time instant t±, the time evolution for the relevant part of G R (ti,t 2 ) 
G^(fi,t 2 ) i = 0, can be calculated effectively. (Note that the integral kernel 

modifies the time derivative of the matrix elements G R 0 , G R 10 only, and e.g., usual 
routines for a set of ordinary differential equations work perfectly.) However, the lower 
limit of the integral in equation (17) still means a difficulty. As a possible way, one 
may use real time decomposition and omit the ’’irrelevant” part of the retarded Green’s 
function as in reference 12 , but in our case there is a more efficient approach that 


completely avoids numerical integration in the Keldysh equation. 


The idea is calculating the time evolution of the integral appearing in equation (26) 
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directly. In order to simplify the notation, let us introduce the following column vector: 


*i(f) = >/N(E)d / G R (t,t')e 


~ iEt 'dt! 


(27) 


The linearity of the Dyson equation allows us to write 

t 

i J^(t) = H(t)*(t) + J - t')<H{t')dt'. (28) 

— OO 

It is worth separating the solution corresponding to H° [the Hamiltonian with no 

(29) 


potential term, see equation (13)], that is, assuming 

^ = T'° + T' 1 , 

where 


*?(«) = VmWjd / 


(30) 


Here, G R0 obeys equation (llTh with U = 0 and the initial conditions given by (|19|). This 


means that G R0 is the retarded Green’s function of the entire potential-free discretized 

infinite line that is evaluated in region II, and consequently it is known analytically: 

G R o(t, t/) = i n Q(t — t')J n [2(t — t')d] exp[— 2i(t — t')d ], where Bessel functions of the 

first kind appear again. Combining this and the analytic (but lengthy) expression for 
o 

f J n (t) exp(iu>t)dt, one obtains 

— OO 

tf°(t) = (3i) 

which is just the continuation of the monochromatic plane wave that arrived from the left 
hand side boundary of region II. (Note the disappearance of the normalization factor.) 
As one can check easily, the equation of motion for T r ° is the following: 

t 

(32) 


iJA°(() = H°»°(t) + J S R (t 


Equation (31) implies that all perturbations induced by the space and time- 
dependent potential U(x,t) is encoded into th 1 . By substracting equation (32) from 


(28), we obtain 


t 

id* 1 (() = U(i)4'“(() + H(()4 ,1 («) + J E fi (( - f')* V)*'. 


(33) 


Additionally, ^ 1 (t) = 0 for t < 0, and it remains zero, unless the potential becomes 


finite; the source term in equation (33) is the first one on the right hand side. 












Quantum mechanical scattering on time-dependent potentials: TDNEGF method 10 


The central equations that provide a solution to the scattering problem we outlined 
in Sec. [2] are (29), (31) and (33). They allow for a clear physical interpretation, for 
which it is the simplest to recall equations 0 and ([26| to see 


n(x = ia, t ) = |dq(f)|' 


(34) 


Additionally, equation (28), which governs the time evolution of ^(t), is essentially a 
Schrodinger equation, apart from points 0 and M, where the integral with memory kernel 
mimics the presence of the semi-infinite line segments I and III, respectively. (That is, 
it provides numerically exact transparent boundary conditions.) Moreover, for U = 0, 
\f '(t) = tE ,0 (f), which is the plane wave input that propagates in an unperturbed way. 

In summary, the complex valued, space and time-dependent function \I/(f) can 
simply be identified with the single particle quantum mechanical wave function. 


5. An application: short laser pulse interacting with an electron beam 

Let us consider an electron beam propagating along the x axis. After t — 0, a short, 
linearly polarized laser pulse impinges on the beam from a perpendicular direction. The 
space and time dependence of the laser field is assumed to be 

x£ 0 cos((nof + </?cep) sin 2 (f7r/T) sin 2 (x7r/L) 
if x £ [0, L\ and t £ [0, T], 


E(x, t) = < 


(35) 


0 otherwise, 

where T characterizes the duration of the pulse. For ultrashort laser pulses, T 
corresponds to few optical cycles only and the carrier-envelope phase <^cep (that 
determines the ’’waveform” of the pulse) can play an important role (X3,fl4l. Note 


of the pulse) can play an important role 13, 14] 
that the spatial and temporal envelopes (sin 2 functions) have been chosen to ensure 
smooth on/off switching. 

We assume that dipole approximation is valid, thus the potential term in the 
Hamiltonian can be written as 


U(x,t) = — / x' E(x' ,t)dx', 


(36) 


where we used atomic units (e = 1). This means that U is not local in the sense of the 
previous sections [see equation ([3])], in region III it is constant in space, but oscillates as 
function of time: Um(t ) = U(L,t)- Consequently, we have to modify our results to be 
able to apply them to this problem. Returning for a moment to the continuous notation 
of the space variable x (and omitting boldface symbols that indicated matrix-vector 
objects), a possible solution is to look for the time-dependent wave function in the form 
of 


T(a;, t ) = u(x, t)4>(x, t ) = u(x , t ) [T°(a;, t ) + 0 1 (x, f)] , 


(37) 
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where u(x,t ) = exp[— i f U(x,t')dt']. Substituting back to equation (28), we obtain: 

o 


0 

i—(f)(x,t) = u*(x,t)H°u(x,t)4>(x,t) 


+ u*(x,t ) / T, R (t — t')(j)(t')u(x, t')dt'. 


(38) 


Subtracting equation (32) (that describes the time derivative of a monoenergetic plane 
wave) from the equation above, an equation of motion for 0 1 (x, t) can be deduced: 

i^(p\x,t) = H°{t)(t) l {x,t) + ~ H°)¥\x,t) 


+ u*(x,t ) / S R (t — t > )(j> 1 (t > )u(x, t')dt', 


(39) 


where H°(t) = u(x,t)H°u*(x,t) and the perfectly transparent boundary property has 
been used to simplify the memory kernel. 



1.01 


1 


0.99 



Figure 3. (Color online) Laser induced density fluctuations in a monoenergetic 
electron beam. The wave function before t = 0 is given by equation (JTJ) , and the 
electron density n(x, t) © is shown in the top row as a contour plot. The difference 
between the columns is the carrier-envelope phase of the driving laser field [see equation 
35 , ycEP = 0 for the left column, </?cep = r/2 for the right one. (As a reference, the 


time dependence of the potential U(x,t ) at x = L is shown in the bottom row.) 


As an example, figure [3] shows the fluctuations of the electron density that arise due 
to the excitation by the laser pulse. Since - according to our remark following equation 
- n(x,t ) = 1 corresponds to the undisturbed plane wave ([!]), we have chosen a 
color scheme where this value is denoted by white (i.e., invisible on the paper). As we 
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expect, the disturbance caused by the laser field becomes visible when the potential gets 
significantly different from zero. The fluctuations, on the other hand, propagate away 
from the central region, and this wave-like behavior is still present when t > T, i.e., when 
the laser pulse is over and the potential is globally zero again. This behavior is related 
to the parameters: the duration of the pulse (35) has the same order of magnitude as 
2 tt/E in figure [3] (concretely: ET=2), which means that the system does not follow 
the change of the potential adiabatically. However, the laser pulse does not mean very 
short, ’’delta-excitation” either, since, as we can see in the figure, the details of the 
time dependence of U(x,t ) strongly modify the function n(x,t), the electron density 
fluctuations are qualitatively different for <pcep = 0 and <Pcep = vr/2 : the system 
exhibits CEP dependence. 


6. Summary 

We described an efficient way of solving scattering problems with time-dependent 
potentials that arc localized both in space and time. We have shown that the theory of 
time-dependent nonequlibrium Green’s functions provides numerically exact dynamical 
equations on a spatial grid so that the boundary conditions are also taken into account 
appropriately. As an example, we have shown that short laser pulses (that contain a 
few optical cycles only) interacting with an electron beam can create electron density 
fluctuations the detailed structure of which depends on the parameters of the exciting 
pulse. 
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